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Abstract 

In solutions of star-branched polyelectrolytes, electrostatic interactions between charged arms 
on neighboring stars can compete with intra-star interactions and rotational entropy to induce 
anisotropy in the orientational distribution of arms. For model stars comprising rigid rodlike arms 
with evenly spaced charged monomers interacting via an effective screened-Coulomb (Yukawa) po- 
tential, we explore the influence of arm orientational anisotropy on effective star-star interactions. 
Monte Carlo simulation and density-functional theory are used to compute arm orientational dis- 
tributions and effective pair potentials between weakly charged stars. For comparison, a torque 
balance analysis is performed to obtain the configuration and energy of the ground state, in which 
the torque vanishes on each arm of the two-star system. The degree of anisotropy is found to 
increase with the strength of electrostatic interactions and proximity of the stars. As two stars 
begin to overlap, the forward arms are pushed back by inter-star arm-arm repulsion, but partially 
interdigitate due to rotational entropy. At center-center separations approaching complete overlap, 
the arms relax to an isotropic distribution. For nonoverlapping stars, anisotropy-induced changes 
in intra- and inter-star arm-arm interactions largely cancel and the effective pair interactions are 
then well approximated by a simple Yukawa potential, as predicted by linear response theory for 
a continuum model of isotropic stars [A. R. Denton, Phys. Rev. E 67, 11804 (2003)]. For overlap- 
ping stars, the effective pair interactions in the simple rigid- arm- Yukawa model agree closely with 
simulations of a molecular model that includes flexible arms and explicit counterions [A. Jusufi, 
C. N. Likos, andH. Lowen, Phys. Rev. Lett. 88, 018301 (2002); J. Chem. Phys. 116, 11011 (2002)]. 
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I. INTRODUCTION 



Polyelectrolytes (PEs) are polymers that carry ionizable groups 1 . In solution, small ions 
(counterions) dissociate, creating oppositely charged polyions. Common examples of PEs 
are polyacrylic acid, sulfonated polystyrene, and biopolymers, such as proteins, DNA, and 
starch. Practical applications in the chemical, microelectronics, and pharmaceutical indus- 
tries include water filtration membranes, direct-write technologies^, DNA-protein binding, 
and thin films for controlled drug delivery 3,4 . When dispersed in an aqueous electrolyte, 
or other polar solvent, PEs interact electrostatically with one another and with charged 
surfaces. Bare Coulomb interactions are screened by dissociated microions (counterions and 
salt ions) distributed in and around the polyions. A fundamental understanding of these mi- 
croscopic interactions is essential for predicting and controlling equilibrium and dynamical 
properties of bulk PE solutions. 

Beyond linear chains, PEs can be readily synthesized with more complex architectures, 
such as stars, microgels, block-copolymer micelles, and dendrimers. In this paper, we focus 
on star-branched PEs, comprising linear PE chains (arms) chemically grafted or adsorbed to 
a common microscopic core. In sharp contrast to colloidal microspheres, whose hard cores are 
rigid and impermeable to microions, PE polyions have internal degrees of freedom associated 
with both chain conformations and microion distributions. Such porous macromolecules are 
characterized by ultra-soft pair interactions^. The main issue addressed here is the extent 
of anisotropy in arm orientations induced by electrostatic interactions between neighboring 
stars and the implications for effective interactions between stars in solution. 

Chain conformations in PE stars have been extensively modeled by scaling theory^, 
self-consistent field theory 7 -, Monte Carlo (MC) 8 and molecular dynamics (MD)*^ simula- 
tions. On the experimental side, several studies have analyzed the structure of PE stars and 
spherical PE brushes in solution. Small-angle neutron scatterin g 11 ' 12 ' 13 and dynamic light 
scattering^ 4 - measurements indicate that the arms of highly ionized stars are almost fully 
stretched, a finding supported by simulations 9 ' 10 . Other neutron scattering studies, using 
isotopically labelled sample a 15 ' 16 ' 17 , have directly probed counterion distributions, indicating 
a strong coupling between counterions and PE chains. These studies, together with osmotic 
pressure measurements^, support conclusions from theoretical and simulation wor k 6 ' 9 ' 10 ' 19 ' 20 
that a high fraction of counterions remain trapped inside the stars, screening electrostatic 
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interactions between the arms. 

Comparatively few studies have examined interactions between PE stars. Recently, we 
applied response theory to a continuum mode l 20 ' 21 , in which the density of charged monomers 
is approximated as an isotropic, continuously varying 1/r 2 distribution, where r is the radial 
distance from the star center. Assuming linear screening of bare Coulomb interactions by 
microions, this approach predicts screened- Coulomb (Yukawa) effective interactions between 
pairs of nonoverlapping stars. In complementary studies, Jusufi, Likos, and Lowen^ - used 
simulation and variational theory to calculate effective pair interactions between overlapping 
stars in a molecular model that includes flexible PE chains and explicit counterions. 

Motivated by recent simulations*^ - and experimen ta 11 ' 12 ' 13 ' 14 , we explore here a relatively 
simple model of PE stars consisting of rigid rodlike arms, lacking bend and twist flexibility, 
and carrying charged beads that interact via screened-Coulomb pair potentials. The intrinsic 
rigidity of PE chains is enhanced in star architecture by electrostatic and solvent-mediated 
excluded-volume repulsions between arms. For example, DNA and some synthetic stiff 
PE chains can be assembled into stars with relatively inflexible arms^ 2 ^ 4 * 2 ^ The rodlike 
model, while directly relevant to such systems, also represents a reference model for stars 
with semiflexible arms. 

The outline of the remainder of the paper is as follows. Section |H] first explicitly defines 
the rigid-arm model of PE stars. Section IIHI then describes three complementary methods 
for analyzing arm orientational anisotropy and effective pair potentials between stars: Monte 
Carlo simulation, torque balance analysis, and density-functional theory. Section HYI next 
presents and discusses numerical results for arm orientational distributions and effective 
pair potentials. Most significantly, the three independent methods yield consistent results 
for star-star interactions, which closely agree with molecular dynamics simulations of the 
molecular model. Finally, Sec. IS summarizes and concludes. 

II. MODEL 

We consider star-branched polyelectrolytes, each star comprising / chains (arms) attached 
to a central core and freely rotating about the core. As depicted in Fig. ^ the arms are 
modeled as rigid rods of equal length a, each arm carrying Nj, charged monomers (beads) 
evenly spaced at a neighbor separation b along the arm. The microions in our model play two 
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essential roles: charge renormalization and electrostatic screening. Most of the counterions 
condense onto or strongly associate with the arms. These bound counterions reduce the 
bare charge on the arms, significantly lowering the effective charge. The remaining mobile 
microions, inside and outside of the stars, screen the bare Coulomb electrostatic interactions. 
Assuming an average effective valence z per bead, the effective valence per star is Z = fN^z. 
The stars are dispersed in an electrolyte solvent, approximated as a dielectric continuum of 
spatially uniform dielectric constant e (primitive model of electrolytes). 

In this paper, we restrict attention to "weakly" charged stars, defined as stars whose arms 
can freely rotate relative to one another and are not orient at ionally localized. As a simple 
physical criterion, a star is weakly charged if the increase in electrostatic energy upon rotat- 
ing an arm half-way towards any nearest neighbor does not significantly exceed the typical 
thermal energy k^T at temperature T. In stars not satisfying this criterion, electrostatic 
interactions overwhelm Brownian motion and the arms are orient ationally frozen. 

To model pair interactions, we focus on two isolated stars whose centers are separated by 
a fixed displacement R. Integrating out the microion degrees of freedo m 2 *- 1 ! 21 and assuming 
linear microion screening of bare Coulomb interactions, pairs of charged beads interact via 
an effective screened-Coulomb (Yukawa) potential. The electrostatic interaction between 
two arms, one oriented along a unit vector u in a star centered at the origin, the other 
oriented along u' in a star centered at R, is then approximated by an effective pair potential 



where e is the proton charge, the summation goes over all the charged beads of the two arms, 
and k is the Debye screening constant (inverse Debye screening length). For a solution with 
star number density p, salt ion pair number density p s , and salt valence z s , the screening 
constant is given by k = y^ATce 2 (z 2 fN b p + 2z 2 p s )/(ek B T). 

The orientational distribution of the arms in two interacting stars is determined by a 
balance of three competing factors. Repulsive forces between arms on different stars exert 
torques that rotate the arms back, creating anisotropy in the arm distribution. Opposing 
this backward rotation, repulsions between arms on the same star exert torques that spread 
out the arms, favoring isotropy. Likewise, random thermal (Brownian) rotational motion 
acts to even out the arm distribution. To calculate the equilibrium arm distribution, and 



z 2 e 2 Nb Nb 
tw(u, u'; R) = z2z2 



exp(— k|R + ibu — jbu'\) 



R + ibu — jbu' 



(1) 
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the corresponding star-star interaction energy, we use three independent methods, described 
in the next section. 

III. METHODS 

A. Monte Carlo Simulation 

For our simple rigid-arm model, Monte Carlo simulation proves to be an efficient method 
for determining the equilibrium distribution of arm orientations at constant temperature. 
The orientation of an arm, denoted by unit vector u, is specified in spherical coordinates by 
polar and azimuthal angles (6, <fi). As illustrated in Fig. |21 for a randomly chosen arm i, a 
trial rotation is performed from the old orientation to a new orientation u\ n ' . Adopting 
an efficient and reliable algorithm^!, we generate a new (trial) orientation according to 

1 11} +7v| 

where v is a unit vector with random orientation and 7 is a tolerance factor that determines 
the magnitude of the trial rotation. The standard Metropolis scheme gives the probability 
of accepting a rotation: 

acc ( ^ n )= min |i ; exp (-/?[$ arm (u( n) ; R) - $ arm (u 4 (o) ; R)]j } , (3) 

where (3 = l/k^T and $ arm (uj;i?) represents the electrostatic energy of the i th arm with 
orientation u« in a star at center-center separation R from a second star. The latter energy 
can be expressed as a sum over arm-arm interactions: 

/ / 
$ arm (u i ; J R)= ^ v axm (u i ,u j ;0) + ^2v aim (u i ,VL j ;R), (4) 

3=1/1 3=1 

where the first summation includes arms of the same star, excluding self-interaction, and 
the second summation includes arms of the second star. 

Our simulations are initialized by centering one star at the origin, the other at a distance 
R away along the z-axis, and assigning random orientations to the arms. Next, trial rotations 
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of the arms are performed in cycles, a cycle consisting of one attempted rotation of each 
arm, selected sequentially. Following an initial equilibration stage of typically 10 3 cycles, 
during which the energy is monitored until it reaches a plateau, arm orientation statistics 
are accumulated in the following 2 x 10 5 cycles. The average energy of the two-star system 
is calculated as 

1 2/ 

$ eff (i?) = -^($ arm ( Uj ;i?)>, (5) 

8=1 

where the summation over all arms includes intra-star and inter-star interaction energies 
and ( ) denotes an ensemble average over configurations. Finally, the effective pair potential 
v e ff(R) between the two stars is obtained as the change in the total energy of the system 
when the two stars are brought from infinite separation to separation R: 

v cS (R) = $eff(i2) - $ cff (oo). (6) 
B. Torque Balance Analysis 

As a check on our simulations, we also calculate the ground-state (T = 0) orientational 
distribution of arms, and the corresponding effective pair potential. At zero temperature, 
where Brownian motion is absent, mechanical equilibrium is reached when the net torque 
on each arm is zero. Electrostatic forces between arms generate torques that drive the arms 
to rotate against the friction of the solvent. Since the terminal angular velocity of an arm 
is proportional to the torque exerted on that arm, the equation of motion is^ 

I&rm 7^ = 7armUj X V Ui $arm(Ui) , (7) 

where J arm and 7 arm are the arm's moment of inertia and rotational friction coefficient, 
respectively. The torque balance analysis amounts to finding the equilibrium configuration 
of the arms such that the right side of Eq. (JZJ) vanishes. 

In practice, we initialize the calculation by assigning random orientations to the arms 
and then use a simple forward- difference scheme to evolve the system until all of the arms 
stop rotating. This procedure efficiently generates the same final distribution as should be 
reached in a Monte Carlo simulation in which the only trial moves accepted are those that 
lower the energy. From the final arm distribution, we compute the ground-state effective pair 
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potential [from Eq. (JHJ)], which represents a lower bound for the pair potential at nonzero 
temperature. 



C. Density- Functional Theory 

An alternative to modeling the arms explicitly is to consider an orientational distribution 
function (ODF) that describes the average configuration of the arms. We define the ODF 
P(u; r) such that P(u; r)du represents the average number of arms with orientation u in 
a solid angle element du subtended at the center of a star at position r. The ODF is 
normalized to the number of arms: J duP(u; r) = /. To approximate the ODF, we propose 
here a simple density-functional theory^ based on a mean-field approximation. 

In the case that translational and rotational time scales are sufficiently separated that 
center-of-mass and orientational coordinates decouple, the grand potential Q[p(r), P(u; r)] 
is a functional of the single-star center-of-mass number density p(r) and the arm ODF 
individually. An arbitrary external potential e xt( r ) uniquely determines p(r) and P(u;r), 
and thus the grand potential^, which can be expressed as 

fi[p(r),P(u;r)] = P id [p(r),P(u;r)] + P cx [p(r),P(u;r)] + J dr J dup(r)P(u; r)[0 ext (r) -p], 

(8) 

where p is the chemical potential of the arms and and P e x are the ideal-gas and excess 
Helmholtz free energy functionals, respectively. The ideal-gas free energy, associated with 
rotational entropy of the arms, is given exactly by 

P id [p(r),P(u;r)] = k B T J dr J dup(r)P(u; r)[ln(47rP(u; r)) - 1], (9) 

translational entropy being ignored, assuming the star centers are fixed. The excess free 
energy, associated with interactions between arms, can be formally expressed as^ 

P ex [p(r),P(u;r)] = ^ J dr J dr' j du J du' J d\ p (2) (r, r'; [Xv arm \) 

x p( 2 )(u,u';r ) r / ;[A^ m ]K m (u,u , ;r-r , ) ) (10) 

where p^(r, r'; [f arm ]) is the two-star density, proportional to the probability of finding at 
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positions r and r' two stars whose arms interact via the pair potential t> arm (u, u';r — r') 
[Eq. ((H)] ; P( 2 )(u, u'; r, r'; [f arm ]) is the two-arm orientational density, proportional to the 
probability of finding two arms at orientations u and u' on stars centered at r and r'; and 
A is a coupling constant that "turns on" the charge. 

For two stars whose centers are fixed and separated by a displacement R, the single-star 
density can be written as 

p(r) = S(r) + 5(r - R) (11) 

and the two-star density as 

p( 2 )(r,r';K m ]) = p(r)p(r'). (12) 

Combining Eqs. (f8J)-()12|). and setting the external potential to zero, the grand potential 
simplifies to 

fi[P(u)] = 2 J P(u){fc B T[m(47rP( U ))-l]-p} 

+ J du J du'^ dA {P^^u'sOjEAVannDVarm 

(u, u'; 0) 

+ F( 2 >(u,u';fi;M)D„(u,u';fi)} 1 (13) 

where the factor of 2 in the first term accounts for the two stars and the dependence of P(u) 
on P is now implicit. In general, the two-arm density can be written as 

P (2) (U, U'j P) = P( U )P(u / )(?ar m (u, u' 5 P), (14) 



which simply defines the arm-arm pair distribution function 5 arm (u, u'; P), in such a way 
that, given an arm with orientation u on a star centered at the origin, P(u')g arm (u, u'; P)du' 
is the average number of arms with orientation u' in solid angle du' on a star centered at 
distance P. 

To facilitate calculations, we make two approximations. First, we neglect arm-arm pair 
correlations and take g a rm(u, u'; P) = 1 for both inter- and intra-star arm-arm pair distri- 
bution functions. This mean-field approximation is based on the assumption that arm-arm 
interactions localize arms to the extent that correlations are mostly accounted for by the 
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anisotropic single-arm distribution function. Because this approximation neglects the self- 
correlation hole at u ~ u', it does give the wrong normalization for intra-star correlations, 
J du J du'P( 2 )(u, u'; 0) = f 2 , compared with the correct result /(/ — 1). However, for multi- 
arm stars this is only a minor inconsistency. In passing, we note that our approximation 
for <7arm( u ) u '; -R) is analogous to that commonly applied to the pair distribution function 
g(r, r') of crystals, whose particles are localized to the extent that the structure of the two- 
particle density, p™(r, r') = p(r)p(r')g(r, r'), is dominated by that of the inhomogeneous 
one-particle density. As a second approximation, we assume that the arms of overlapping 
stars (R < 2a) do not interdigitate, but instead rotate backwards, being strictly excluded 
from a forward cone (Fig. |HJ), within which P(u) = 0. The validity of the no-interdigitation 
approximation is examined below in Sec. II VI 

With the above two approximations, Eq. (|13|) simplifies to 

n[P(u)] = 2 J duP(u){k B T[\n(AnP(u)) - 1] - fi} 

+ J du J du / P( U )P(u / )Km(u,u , ;0)+t; arm (u,u / ; J R)]. (15) 

At equilibrium, fi[P(u)] is a minimum with respect to P(u). Taking a functional derivative 
with respect to P(u), we obtain 

15fig(u)] = fcBTln(47rP(u)) + f du 'p( u ')[,; arm ( U)U ' ; o) +t, arm (u,u';P)]. (16) 
2 dP{u) J 

Setting the right side of Eq. (fTfiJl to zero yields an implicit equation for P(u): 

P(u) = ^exp (-P J du / P(u')K rm (u,u / ;0) +^ arm (u,u , ;P)]^ . (17) 

On the right side of Eq. ([170. the first term of the integrand corresponds to interactions 
between arms in the same star (intra-star interactions) and the second term to interactions 
between arms in different stars (inter-star interactions). Note that, by symmetry, the two 
stars are mirror images of each other reflected in the perpendicular plane bisecting the line 
joining the star centers. 

Once the equilibrium orientational distribution of arms is determined, the energy of the 
two-star system - a sum of the intra-star and inter-star arm interaction energies - can be 



9 



obtained from 



$ eS {R) = J du J du / P( U )P(u / )Km(u,u / ;0)+t; arm (u,u , ; J R)], (18) 

which is the continuum analog of Eqs. and (0) for the discrete model. The effective 
star-star pair interaction is finally calculated via Eq. (jSj). taking P — > f/An as i? — > oo. 



IV. RESULTS AND DISCUSSION 

Having described three methods for calculating the arm orientational distribution and the 
effective star-star pair potential, we now present and discuss numerical results. Our choice 
of parameters is restricted here by the defining criterion for weakly charged stars (Sec. IrT]) . 
In this parameter regime, where azimuthal arm-arm correlations are weak, we can exploit 
the axial symmetry of the ODF with respect to the line joining the centers of the two stars. 
Thus, the ODF is treated as a function of only the polar angle: P(u) = P{0). We have 
checked and confirmed this assumed symmetry in our MC simulations. 

To extract P(0) from the simulations, we bin the arm orientations in a manner illustrated 
in Fig. I3J Each sphere is partitioned into slices perpendicular to the line joining the star 
centers. Allowing each slice to subtend the same angle A9, the average number of arms 
in the k th slice centered at angle 9k is 2nP(9k) sin 6k AO. The same number can also be 
expressed as fN^J Ylh=i wnere is the number of occurrences (accepted trial moves) 
of an arm in the k th slice and M is the total number of slices. Equating the two expressions, 

P(e k ) = fNk M — , (19) 

27rsin0 fe A0£f =1 ivV 
which is correctly normalized to the number of arms /. 

To compute the corresponding DFT approximation for P(0), we solve Eq. (|17|) iteratively, 
computing the arm-arm pair potentials t> arm (u, u'; R) from Eq. (JTJ. In practice, we construct 
a 60 x 60 grid on a sphere, with 60 equi-spaced longitudinal lines in the azimuthal coordinate 
(0 < cf) < 2tt) and 60 equi-spaced latitudinal lines in the polar coordinate (0 < < 7r). 
Starting from an initial normalized isotropic distribution, P = //47r, we solve Eq. ()17)) self- 
consistently by numerical iteration until convergence. As a numerical check, we obtain the 
same results using Monte Carlo integration to compute the angular integral in Eq. ()17|) by 
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summing over randomly sampled directions. 

Figure 03 compares our MC simulation results and DFT predictions for 8-arm stars, each 
of radius a = 10 nm and with 10 beads per arm, for a Debye screening constant Ka = 1 
and various star valences and separations. The error bars represent root-mean-squared 
deviations from the average of intermediate averages over 100-cycle intervals. Evidently, the 
higher the valence and the closer the separation of the two stars, the more pronounced the 
anisotropy of the ODF, consistent with intuition. In the case of nonoverlapping stars, the 
DFT predictions agree well with simulation for a relatively low valence (Z = 20), as seen 
in Fig. Efa). For a higher valence [Z = 55), the mean- field theory's neglect of arm-arm 
correlations leads to an underestimation of the ODF structure [Fig. E(b)]. 

In the case of overlapping stars, the DFT's simplifying assumption of no interdigitation 
qualitatively models suppression of the ODF in the forward (small-6 1 ) direction, but the 
sharp cut-off in P(9) is clearly too severe. As Fig. |^c) shows, Brownian motion at nonzero 
temperature can cause the forward arms of two overlapping stars to partially interdigitate. 
To quantify this behavior, we define an interdigitating configuration as one in which an arm 
on one star pierces a triangle whose vertices are the tips of two arms on the second star and 
the center of that star. Based on this criterion, we define the "interdigitation ratio" as the 
number of accepted configurations with interdigitating arms divided by the total number of 
accepted configurations. As Fig. El shows, the interdigitation ratio increases with decreasing 
star separation and with decreasing star valence. Nevertheless, the extent of interdigitation 
is sufficiently weak that the DFT's neglect of such configurations still yields an accurate 
effective pair potential (see below). 

To further quantify the extent of anisotropy in the arm orientational distribution, we 
define an orientational order parameter 



For a perfectly isotropic distribution, S — 0, whereas for a highly anisotropic distribution for 
which the ODF peaks at 9 = it, S = 1. Extracting the order parameter from a simulation 




(20) 
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simply amounts to computing the sum 

1 - 

S=-— B —Y,*k<xxO u . (21) 

2^=1 fc=1 

Figure [7| displays the orientational order parameter as calculated from simulation for vari- 
ous choices of parameters, namely star valence, arm number, and screening length. As seen 
in Figs. [7^ a) and (b), for fixed arm number the order parameter increases with increasing 
star valence and increasing screening length (decreasing screening constant). Figure [TJc) 
demonstrates that for fixed valence and screening length the order parameter increases 
upon reducing the number of arms (increasing charge per arm). Thus, orientational or- 
dering (anisotropy) is amplified by strengthening or increasing the range of electrostatic 
interactions, which is consistent with the trends observed in the ODF (Fig. EJ). 

Interestingly, the degree of anisotropy does not vary monotonically with center-center 
distance between two stars. On the contrary, it attains a maximum at a distance where two 
stars strongly overlap, but short of complete overlap. The abrupt decrease in S at very short 
separations simply reflects the loss of identity of completely overlapping stars that merge into 
one isotropic "super-star." As an aid to visualizing the anisotropy, Fig. |S] provides snapshots 
of the arm orientations of two 8-arm stars, as determined from the torque balance analysis 
and MC simulation. 

Next, we compare results for the effective star-star pair potential, as computed from MC 
simulation, torque balance analysis, and DFT [Eqs. (fTTj) . and (|TK|)]. As seen in 

Fig. El the torque balance analysis predicts the lowest pair energy of the three methods, as 
expected of a method designed to find the ground state configuration. The actual energy 
is higher because of thermal rotational motion of the arms at nonzero temperature. The 
DFT predictions closely track the simulation data, aside from relatively small deviations at 
overlapping separations. 

The consistent agreement between theory and simulation at all separations is somewhat 
unexpected, considering the significant discrepancies in arm orientational distributions for 
high valences and for overlapping stars [Fig. Efc)]. The theory's complete neglect of arm 
correlations and interdigitation is obviously unphysical for significantly overlapping stars, 
especially in the forward direction, i.e., for small 9. As a consequence, the theory underesti- 
mates the density of (interdigitating) arms within the forward cone [Fig. [2] and overestimates 
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the arm density just outside of the forward cone. For stars that are only weakly overlapping, 
the ODF discrepancies are limited to arm configurations with small polar angles, which ac- 
count for a relatively low fraction of the phase space in the angular integral for the effective 
pair potential [Eq. (|18|) ]. For more strongly overlapping stars, the close agreement between 
theory and simulation must be attributed to a fortuitous cancellation of errors - namely 
an overestimate of the contribution to v c r (R) from arms outside of the interdigitation cone 
compensated by an underestimate of the contribution from interdigitating arms. Although 
it should be feasible to devise an analytical fit to our results for the effective pair potential 
of two overlapping stars for any set of system parameters, namely star radius, arm number, 
linear charge density, and screening length, we leave this fitting analysis to a future study. 

For nonoverlapping stars, the simulation and DFT results hardly deviate from the pre- 
dictions of linear response theory applied to an isotropic continuum model for the 1/r 2 
charge distribution of PE stars 2 ^. The latter theory predicts a screened- Coulomb (Yukawa) 
potential between pairs of stars of the form 



v e s{R) 



Z 2 e 



2^2 



shi(/ta) ( 



KCL 



2 a - K R 



R 



R > 2a, (22) 



where shi(x) = f Q du smh(u)/u denotes the hyperbolic sine integral function. The accuracy 
of Eq. (}2"2*|) may seem surprising, given the anisotropic arm distributions (Fig. EJ). Closer 
examination reveals, however, that the anisotropy-induced reduction of inter-star interaction 
energy is largely balanced by an increase of the intra-star interaction energy, as illustrated 
in Fig. PP . 

Finally, as a test of the validity of the rigid- arm- Yukawa model and the accuracy of our 
methods, we compare with results from molecular dynamics simulations of a molecular model 
of PE stars&i£ in which the arms are represented as bead-spring chains of Lennard- Jones 
particles connected by a finite-extension-nonlinear-elastic (FENE) potential and charged 
monomers and counterions interact via bare Coulomb potentials. The strong stretching of 
chains observed in the MD simulations^^ allows direct comparison with our simpler rigid- 
arm model. To facilitate comparison, we choose parameters to match those given in Table 
IV of Jusufi et a/.—, including the Debye screening length, Bjerrum length e 2 jek^T = 0.75 
nm, and fraction of charged beads per arm a. Following refs.^ 10 , we treat the effective star 
valence Z as a tunable parameter to fit the effective pair force F = —dv e g(R) / OR. (Note that 
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tuning Z also changes the screening length.) From the effective valence, we then obtain the 
effective number of counterions condensed inside the two stars, N\ = 2(fNi,a — Z). As shown 
in Fig. El our best fits to the MD forces between overlapping pairs of stars with charging 
fraction a = 1/4 are obtained by choosing Ni = 137 for 10-arm stars and Ni = 465 for 30- 
arm stars. For comparison, the corresponding MD values^ are, respectively, iV\ = 147 and 
Ni = 450. The excellent fit of the MC and MD force data, in particular the close agreement 
of the effective number of condensed counterions, provides a consistency check and supports 
the simple rigid- arm- Yukawa model for the weakly charged stars considered here. Further 
such comparisons would be useful to establish the limits of validity of the model. 

Three issues merit further discussion. First, while molecular models that include chain 
flexibility and explicit microions naturally provide the most realistic description of PE stars, 
the many degrees of freedom present major computational challenges. Simulations of such 
models have thus far been limited to only one or two star a 8 i 9 i 10 . By comparison, the rela- 
tively coarse-grained model explored here seems to accurately capture effective interactions 
between PE stars, yet is simple enough to allow studies of phase behavior in simulations of 
many interacting stars. 

Second, although simulations of the coarse-grained and molecular models yield very sim- 
ilar pair forces, the rather different treatments of the microions in the two models suggest 
complementary physical interpretations. In the molecular model, the interactions between 
overlapping stars are attributed largely to entropy of trapped counterions, electrostatics 
playing only a relatively minor rol e 9 ' 1 ' - ' . In contrast, in the coarse-grained one-component 
model, where the microions are integrated out in a first step to obtain effective (screened- 
Coulomb) arm-arm interactions, electrostatics is the sole determiner of effective star-star 
interactions. Thus, the entropy of trapped counterions, which in the molecular model is 
explicitly associated with the distribution of counterions around the arms of the stars, has 
its coarse-grained counterpart in the implicit contribution of counterions to the effective 
star valence and to the effective interactions via the Debye screening constant. Although 
the correspondence between the two models has been tested here using the available MD 
data for 10-arm and 30-arm stars [Fig. Ill] , it should be tested more extensively in future, 
especially over a broader range of star functionality. 

Third, the PE stars considered here belong to a class of macromolecules for which chain 
entropy is negligible. Since we model the arms as charged rigid rods, and ignore excluded- 
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volume interactions between arms, the effective pair interactions vanish in the limit of charge- 
neutral stars. For stars with flexible arms, Witten and Pincus^ and Likos et alM have shown 
that conformational entropy of the arms generates an effective "ultrasoft" repulsion between 
neutral stars, which varies logarithmically with distance at small separations and follows a 
Yukawa form at longer separations. Such entropic interactions are here neglected. 

V. SUMMARY AND CONCLUSIONS 

In summary, within a coarse-grained model of polyelectrolyte stars, we have analyzed arm 
orientational anisotropy and effective interactions between an isolated pair of weakly charged 
(orient at ionally molten) stars. The arms of the stars were modeled as rigid rods, freely ro- 
tating about the star centers, and carrying evenly spaced charged beads. Microions were 
included implicitly by modeling electrostatic interactions between pairs of charged beads 
via an effective screened- Coulomb (Yukawa) potential. Arm orientational distribution func- 
tions and effective star-star electrostatic interactions were calculated by three independent 
methods: Monte Carlo simulation, torque balance analysis, and density-functional theory. 

The arm orientational distributions predicted by DFT are in good agreement with sim- 
ulation for nonoverlapping stars of relatively low effective valence. With increasing overlap 
and increasing valence, the agreement worsens due to the theory's neglect of arm-arm cor- 
relations and interdigitation. An orientational order parameter, quantifying the extent of 
anisotropy, increases as two stars approach each other, peaks at strong overlap, and then 
decreases as the stars near complete overlap. 

The simulation and DFT results for the effective pair potential between nonoverlapping 
stars are in excellent agreement with previous predictions of linear response theory applied 
to a continuum model of PE stars with isotropic 1/r 2 charge distribution^. This close 
agreement results not only from the relatively weak anisotropy of separated stars, but also 
from a fortuitous cancellation between reduced inter-star interactions and enhanced intra- 
star interactions induced by arm anisotropy. We conclude that electrostatic interactions 
between nonoverlapping weakly charged PE stars can be accurately modeled by a simple 
Yukawa effective pair potential. 

The effective pair potentials from our Monte Carlo simulations and DFT calculations are 
in good agreement over the whole range of star separations, including strongly overlapping 
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stars, for the parameters here investigated. Furthermore, treating the effective valence of 
the st fitting parameter, we achieve a close fit to available data for the effective pair 

forces and numbers of condensed counterions from MD simulations*^ of a molecular model 
with flexible bead-spring chains and explicit counterions. 

Simulations of many-star systems are a natural next step. For dilute solutions, inter- 
actions between isolated pairs of nonoverlapping stars can be reliably approximated, as 
confirmed here, by an isotropic screened- Coulomb pair potential 2 ^. Only interactions be- 
tween overlapping stars must then be computed by explicit summation over pairs of charged 
beads, pending an analytical fit to our MC data. For concentrated solutions of frequently 
overlapping stars, arm anisotropy and many-body interactions among three or more stars 
are likely to be important. In this regime, explicit pairwise summation over beads is es- 
sential. Such an approach could independently test the phase behavior predicted by recent 
simulations of dense solutions of stars interacting via isotropic effective pair potentials 34 . 
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(a) (b) 

FIG. 1: (a) Schematic drawing of PE star and counterions. (b) Model of PE star with rigid, 
rodlike arms. Integrating out degrees of freedom of the microions leads to an effective Yukawa pair 
potential between pairs of charged beads along the arms. 




FIG. 2: Trial rotation of the i th arm from its old orientation u- to its new orientation u- . 
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FIG. 3: Geometry assumed in density-functional theory of overlapping stars. Arms are completely 
excluded from cone-shaped volumes (white region). Darker shading indicates higher arm density. 




FIG. 4: Illustration of Monte Carlo data analysis. Dashed lines represent slices of spheres perpen- 
dicular to center-center line between stars. Points indicate positions of arm tips at various polar 
angles 9, each slice subtending an angle A9. Statistical analysis of configurations yields the arm 
orientational distribution function. 
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FIG. 5: Arm orientational distribution function P(9) vs. polar angle 6 for two 8-arm stars of radius 
(arm length) a = 10 nm, with 10 beads per arm, Debye screening constant kcl = 1, and various 
labeled star valences Z and center-center separations R. Results from Monte Carlo simulation 
(symbols) and density-functional theory (curves) are compared. Error bars represent statistical 
uncertainty (one standard deviation) in simulation data. 

20 




R/a 



FIG. 6: Fraction of interdigitating configurations (see text) for two overlapping 8-arm stars of 
various valences at fixed screening constant na = 1. 
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FIG. 7: Orientational order parameter S vs. center-center separation R between two stars of radius 
a = 10 nm, with 10 beads per arm, for various parameter combinations, (a) Fixed arm number 
and star valence, varying screening length, (b) Fixed arm number and screening length, varying 
star valence, (c) Fixed screening length and valence, varying arm number. Symbols are simulation 
data and curves are guides to the eye. 



FIG. 8: Configurations of two approaching stars at center-center distance R/a = 1 (upper panel) 
and R/a = 0.2 (lower panel) from a torque balance analysis. Each star has 8 arms of length 10 
nm, with 10 beads per arm, valence Z = 20, and screening constant ko, = 1. 
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FIG. 9: Effective pair potential vs. center-center separation between an isolated pair of 8-arm 
stars of valence (a) Z = 20 and (b) Z = 55, calculated from Monte Carlo simulation (•), density- 
functional theory (solid line), torque balance analysis (A), and linear response theory+continuum 
model 20 (dashed line). Insets magnify nonover lapping region. 
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FIG. 10: Monte Carlo simulation results for anisotropy-induced changes in inter- and intra-star pair 
interaction energies for 8-arm stars with arm length a = 10 nm, 10 beads per arm, and screening 
constant na = 1. Solid and dashed curves are intra-star (upper) and inter-star (lower) energies for 
star valences Z = 20 and Z = 60, respectively. Inset magnifies nonoverlapping region. 
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FIG. 11: Comparison of forces between two stars of radius a vs. center-center separation R 
from our Monte Carlo (MC) simulations of the rigid-arm- Yukawa model (curves) and molecular 
dynamics (MD) simulations of a molecular model (symbols) 10 . Choosing the number of condensed 
counterions as N\ = 137 for 10-arm stars and N\ = 465 for 30-arm stars yields good agreement 
between MC and MD data (see text). 
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